//============================================================================
// Name        : FabVoronoi.cpp
// Author      : 
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <iterator>
#include <iomanip>
#include "Topol.h"
#include "MAtoms.h"
#include "HeaderTrj.h"
#include "ResidueCM.h"
#include "FComms.h"
#include "Help.hpp"
#include "FstreamC.h"
#include "FstreamF.h"
#include "trjInput.h"
#include "Gyration.h"
#include "RhoStat.h"

using namespace Topol_NS;
using namespace std;
unsigned int nnx=128,nny=128,nnz=128;
// profile variables
string sVect="X";
double ddx1=0.1*unit_nm; // in Angstroems.
double ddx2=1.5*unit_nm; // in Angstroems.
int psave=20;
double cutX=5.0*unit_nm;
double cutYZ=40.0*unit_nm;
bool bprof=false;

int main(int argc, char ** argv)
{
	trj::trjInput MyInput(argc,argv);
	Parallel::comm=new Parallel::MyMPI(argc,argv);
	bool bDel=true;
	bool bWriteFab=false;
	bool bTestVol=false;
    bool bHyd=true;
	int narg=1;
	string filein,fileout="Micelles.out",filepdb,filedomain,filedefdomain,fileirl,fileidb;
	string fileoutp1,fileoutp2,fileoutp3;

	string scratchdir;
	ifstream fpdb;
	ifstream ftest,fdefdomain;
	Fstream * finx=NULL;
	Fstream * fout_xtcx=NULL;
	ofstream * fout_pdbx=NULL;
	ofstream * foutx=NULL;
	ofstream * foutp1=NULL;
	ofstream * foutp2=NULL;
	ofstream * foutp3=NULL;
	ofstream fdomain;
	string comm0(argv[0]);
	size_t mypos=comm0.find_last_of("/")+1;
	size_t length=comm0.size()-mypos;
	string command=comm0.substr(mypos,length);
	string errmsg;
	string Usage="Usage:\t"+ command + "\n";
	vector<string> use=MyInput.getUsage();
	vector<string> SelRes;
	string Reference;
	for(unsigned int n=0;n<use.size();n++)
		Usage+=use[n];
	Usage+="\n\t Default values in square brackets []\n";
	try{
		if(argc == 1) throw Usage;
		if(int m=MyInput.bTest().size()) {
			errmsg=" Command(s) not found: ";
			for(unsigned int n=0;n<m;n++)
				errmsg+=MyInput.bTest()[n]+"  ";
			errmsg+="\n"+Usage;
			throw errmsg;
			}
		}
	catch(const string & s){
		cout << s << endl;
		Parallel::comm->Barrier();
		Parallel::comm->~MyMPI();
		return 0;
	}


	int nstart=0;
	int nend=-1;
	int nskip=1;
	// input read
	bool inputfile=false;
	try{
		if(!MyInput["-dcd"].empty()) {
			if(MyInput["-dcd"].size() < 2) throw string("\n Filename expected for " + MyInput["-dcd"][0] + " option \n");
			if(MyInput["-dcd"].size() > 2) throw string("\n More than one entry for " + MyInput["-dcd"][0] + " option \n");
			filein=MyInput["-dcd"][1];
			ftest.open(filein.c_str(),ios::in);
			if(!ftest) throw string("\n Cannot open " + filein + "!!\n");
			ftest.close();
			finx=new FstreamF(filein);
			inputfile=true;
		}
		if(!MyInput["-xtc"].empty()) {
			if(MyInput["-xtc"].size() < 2) throw string("\n filename expected for " + MyInput["-xtc"][0] + " option\n ");
			if(MyInput["-xtc"].size() > 2) throw string("\n More than one entry for " + MyInput["-xtc"][0] + " option \n");
			filein=MyInput["-xtc"][1];
			ftest.open(filein.c_str(),ios::in);
			if(!ftest) throw string("\n Cannot open " + filein + "!!\n");
			ftest.close();
			finx=new FstreamC(filein,"rb");
			inputfile=true;
		}
		if(!MyInput["-pdb"].empty()) {
			if(MyInput["-pdb"].size() < 2) throw string("\n filename expected for " + MyInput["-pdb"][0] + " option \n");
			if(MyInput["-pdb"].size() > 2) throw string("\n More than one entry for " + MyInput["-pdb"][0] + " option \n");
			filepdb=MyInput["-pdb"][1];
			fpdb.open(filepdb.c_str(),ios::in);
			if(!fpdb) throw string("\n Cannot open " + filepdb + "!!\n");
		}
		if(!MyInput["-select"].empty()) {
			if(MyInput["-select"].size() != 2) throw string(" String of selected residues needed for " + MyInput["-select"][0] + " option ");
			string selection=MyInput["-select"][1];
			stringstream iss(selection);
			copy(istream_iterator<string>(iss),
					istream_iterator<string>(),
					back_inserter<vector<string> >(SelRes));
		}
		if(!MyInput["-center"].empty()) {
			if(MyInput["-center"].size() != 2) throw string(" Reference residue for Micelle centering needed " + MyInput["-center"][0] + " option ");
			Reference=MyInput["-center"][1];

		}
		if(!MyInput["-o"].empty()) {
			if(MyInput["-o"].size() < 2) throw string("\n filename expected for " + MyInput["-o"][0] + " option \n");
			if(MyInput["-o"].size() > 2) throw string("\n More than one entry for " + MyInput["-o"][0] + " option \n");
			fileout=MyInput["-o"][1];
		}
		if(!MyInput["-b"].empty()) {
			if(MyInput["-b"].size() != 2) throw string(" Number of first frame needed for " + MyInput["-b"][0] + " option ");
			stringstream(MyInput["-b"][1])>> nstart;
		}
		if(!MyInput["-e"].empty()) {
			if(MyInput["-e"].size() != 2) throw string(" Number of last frame needed for " + MyInput["-e"][0] + " option ");
			stringstream(MyInput["-e"][1])>> nend;
		}
		if(!MyInput["-skip"].empty()) {
			if(MyInput["-skip"].size() != 2) throw string(" Number of skipped frames needed for " + MyInput["-skip"][0] + " option ");
			stringstream(MyInput["-skip"][1])>> nskip;
		}
		if(!MyInput["-nohyd"].empty()) {
			if(MyInput["-nohyd"].size() != 1) throw string(" No argument for " + MyInput["-nohyd"][0] + " option ");
				bHyd=false;
		}
		if(!MyInput["-hyd"].empty()) {
			if(MyInput["-hyd"].size() != 1) throw string(" No argument for " + MyInput["-hyd"][0] + " option ");
				bHyd=true;
		}
		if(!MyInput["-nodel"].empty()) {
			if(MyInput["-nodel"].size() != 1) throw string(" No argument for " + MyInput["-nodel"][0] + " option ");
				bDel=false;
		}
		if(!MyInput["-del"].empty()) {
			if(MyInput["-del"].size() != 1) throw string(" No argument for " + MyInput["-del"][0] + " option ");
				bDel=false;
		}
		if(!MyInput["-prof"].empty()) {
			if(MyInput["-prof"].size() != 1) throw string(" No argument for  " + MyInput["-prof"][0] + " option ");
			bprof=true;
		}

		if(!MyInput["-cut"].empty()) {
			if(MyInput["-cut"].size() != 2) throw string(" General cutoff needed for  " + MyInput["-cut"][0] + " option ");
			stringstream(MyInput["-cut"][1])>> cutYZ;
			cutYZ*=unit_nm;
		}
		if(!MyInput["-cutX"].empty()) {
			if(MyInput["-cutX"].size() != 2) throw string(" Out of plane cutoff needed for  " + MyInput["-cutX"][0] + " option ");
			stringstream(MyInput["-cutX"][1])>> cutX;
			cutX*=unit_nm;
		}
		if(!MyInput["-pdx1"].empty()) {
			if(MyInput["-pdx1"].size() != 2) throw string(" Bin size needed for  " + MyInput["-pdx1"][0] + " option ");
			stringstream(MyInput["-pdx1"][1])>> ddx1;
			ddx1*=unit_nm;
		}
		if(!MyInput["-pdx2"].empty()) {
			if(MyInput["-pdx2"].size() != 2) throw string(" Bin size needed for  " + MyInput["-pdx2"][0] + " option ");
			stringstream(MyInput["-pdx2"][1])>> ddx2;
			ddx2*=unit_nm;
		}
		if(!MyInput["-psave"].empty()) {
			if(MyInput["-psave"].size() != 2) throw string(" Bin size needed for  " + MyInput["-psave"][0] + " option ");
			stringstream(MyInput["-psave"][1])>> psave;
		}
		if(!MyInput["-grid"].empty()) {
			if(MyInput["-grid"].size() != 4) throw string(" Grid points in x,y,z needed for  " + MyInput["-grid"][0] + " option ");
			stringstream(MyInput["-grid"][1])>> nnx;
			stringstream(MyInput["-grid"][2])>> nny;
			stringstream(MyInput["-grid"][3])>> nnz;
		}
		if(!MyInput["-pVct"].empty()) {
			if(MyInput["-pVct"].size() != 2) throw string(" Perpendicular axes needed for  " + MyInput["-pVct"][0] + " option ");
			sVect=stringstream(MyInput["-pVct"][1]).str();
		}
		if(!MyInput["-help"].empty()) {
			if(MyInput["-help"].size() != 1) throw string(" No argument for " + MyInput["-help"][0] + " option ");
				throw Help;
		}
	}catch(const string & s){
		cout << s << endl;
		Parallel::comm->Barrier();
		Parallel::comm->~MyMPI();
		return 0;
	}
	if(fileout.substr(fileout.find_first_of(".")+1).find("xtc") != string::npos){
		fout_xtcx=new FstreamC(fileout,"wb");
	} else if(fileout.substr(fileout.find_first_of(".")+1).find("pdb") != string::npos){
		fout_pdbx=new ofstream(fileout.c_str(),ios::out);
	} else if(fileout.substr(fileout.find_first_of(".")+1).find("prf") != string::npos){
		bprof=true;
	} else {
		foutx=new ofstream(fileout.c_str(),ios::out);
	}

	if(bprof){
		string tmp1=fileout.substr(0,fileout.find_first_of("."));
		string tmp2=fileout.substr(fileout.find_first_of("."),fileout.size()-1);
		fileoutp1=tmp1+"-1D"+tmp2;
		fileoutp2=tmp1+"-2D"+tmp2;
		fileoutp3=tmp1+"-3D"+tmp2;
		foutp1=new ofstream(fileoutp1.c_str(),ios::out);
		foutp2=new ofstream(fileoutp2.c_str(),ios::out);
		foutp3=new ofstream(fileoutp3.c_str(),ios::out);
	}
	Topol::fDomainPDB=new std::ofstream ;
	Topol::fDomainPDB=&fdomain;

	RhoStat::setNN(nnx,nny,nnz);
	RhoStat Rho_d;
	sVect.erase(remove_if(sVect.begin(),sVect.end(),::isspace),sVect.end());
	char vect=sVect.c_str()[0];
	Rho_d.sHistAxis(vect);

	Fstream & fin=*finx;
	Fstream & fout_xtc=*fout_xtcx;
	ofstream & fout_pdb=*fout_pdbx;
	ofstream & fout=*foutx;
	ofstream & foutpp1=*foutp1;
	ofstream & foutpp2=*foutp2;
	ofstream & foutpp3=*foutp3;

	try{
		if(!inputfile) throw string("\nNo input file given will operate on the .pdb file\n");
	}
	catch(const string & s ){
		cout << s << endl;
		Parallel::comm->Barrier();
	}

	// read pdb file to construct topology
	vector<string> data;
	for(string str;getline(fpdb,str);){
		data.push_back(str);
	}
	Topol MyTop;

	MyTop(data,false);

	typedef map<string,vector<vector<int> > > mappa;
	mappa & Def=MyTop.getDef();
	try{
		if(Reference.empty()) throw string("Reference residue for centering is needed in argument [-center]. This parameter is mandatory."
				"\n Please input reference residue now: ");
		else{
			if(!Def.count(Reference)) throw string("Unknown reference residue "+Reference+
					"\n Please input reference residue now: ");
		}
	} catch(const string & s){
		cout << "\n System residues are: ";
		for(mappa::iterator it=Def.begin();it != Def.end();it++){
			cout << it->first << " " ;
		}
		cout<< "\n" << endl;
		cout << s ;
		bool ok=true;
		while(ok){
			cin >> Reference;
			if(cin.eof()) {
				cout << " \n ctrl-D character entered. Exit now !\n" <<endl;
				exit(1);
			}
			if(!Def.count(Reference)) {
				cout << "Unknown reference residue "+Reference+
					"\n Please input reference residue now:\n % ";
			}
			else ok=false;
		}
	}

	vector<string> CSelect;
	if(SelRes.empty()){
		cout << "\n System residues are: ";
		for(mappa::iterator it=Def.begin();it != Def.end();it++){
			cout << it->first << " " ;
		}
		cout << endl;
		cout << "\nSelection of the micelle residues not given. \nPick "
				"the residues now  and type ctrl-D to finish: ";
	} else {
		for(unsigned int n=0;n< SelRes.size(); n++)
			if(!Def.count(SelRes[n])) CSelect.push_back(SelRes[n]);
		if(!CSelect.empty()){
			cout << "\n System residues are: ";
			for(mappa::iterator it=Def.begin();it != Def.end();it++){
				cout << it->first << " " ;
			}
			cout << endl;
			cout << "\n Residue(s) " ;
			for(unsigned m=0; m<CSelect.size();m++)
				cout << CSelect[m] << "  ";
			cout <<" are not on the list.\n\nPick "
				"the residues now  and type ctrl-D to finish:\n";
		}

	}
	if(SelRes.empty() || CSelect.size()){
		SelRes.clear();
		bool ok=true;
		do {
			ok=true;
			string dummy;
			bool po=true;
			while(cin>>dummy) {
				SelRes.push_back(dummy);
			};

			for(unsigned int n=0;n < SelRes.size();n++)
				if(!MyTop.CheckResidue(SelRes[n])) {
					ok=false;
					cout << SelRes[n] << " " ;
				}

			if(ok) {
				cout << "\nSelected Residues are:  ";
				for(unsigned int n=0;n < SelRes.size();n++)
				cout << SelRes[n] << " ";
				cout <<endl;
			} else{
				cout << " not on the list. Try again.\n";
				SelRes.clear();
				cin.clear();
			}
		} while(!ok);
	}

// Micelle stuff

		int natoms=MyTop.Size();
	int nresid=MyTop.ResSize();
	MAtoms atm;
	atm.setDim(natoms);

	ios::streampos len;
	HeaderTrj header;
// Read header of dcd file
	if(finx) {
		fin.seekg(0,"end");
		len=fin.tellg();
		fin.seekg(0,"beg");
		fin>>header;
		try{if(!header.check(MyTop.Size())) throw string("Number of atoms in the pdb and trajectory files does not match!");}
		catch(const string & s){cout << s<<endl;exit(-1);}
	}

	atm.initPDB(MyTop.getPDB());

	atm.setList(SelRes);
// Read trajectory and compute Micelle data

    Parallel::comm->Barrier();

    Histogram3D plot3D;
    Histogram2D plot2D;
    Histogram1D plot1D;
    if(!finx){
    	atm.pdb(data);
    	atm.Reconstruct(Reference,MyTop);
      	if(fout_xtcx) fout_xtc << atm ;
    	if(fout_pdbx) fout_pdb << atm ;
    	if(foutx) {
    		atm.Gyro();
    		fout << atm.getRg_i();
    	}
    	Parallel::comm->Barrier();
    	Parallel::comm->~MyMPI();
    	cout << "\nProgram completed: Output data written to " + fileout << "\n\n";
        return 0;
    }
	int ntime=0;
	int nnStep=0;
	while(fin.tellg() < len){
		if(ntime < nstart) {fin+=atm;ntime++;continue;}
		if(nend != -1 && ntime > nend) break;
		fin >> atm;
		cout << fixed << setw(5) << "----> Step No. " << setw(5) << right << fixed << ntime << " -- Time = " << atm.getTime() <<" ps \n";
    	atm.Reconstruct(Reference,MyTop);
    	if(fout_xtcx) fout_xtc << atm ;
    	if(fout_pdbx) fout_pdb << atm ;
    	if(bprof) Rho_d.DensityRot(atm);
    	if(foutx){
    		atm.Gyro();
    		fout << atm.getRg_i();
    	}
    	Parallel::comm->Barrier();
		ntime++;
		nnStep++;
	};
	if(bprof){
		Rho_d.doStat();
		plot3D=Rho_d.Density3D(ddx2,cutYZ);
		plot2D=Rho_d.DensityP2D(ddx2,ddx2,cutYZ,cutX);
		plot1D=Rho_d.DensityP1D(ddx1,cutYZ);
		foutpp3 << plot3D;
		foutpp2 << plot2D;
		foutpp1 << plot1D;
	}

	Parallel::comm->Barrier();
	Parallel::comm->~MyMPI();
	if(bprof){
		cout << "\nProgram completed: Output data written to files: " + fileoutp1+" and "+ fileoutp2  << "\n\n";
	} else cout << "\nProgram completed: Output data written to " + fileout << "\n\n";
    return 0;
}   //end of main function
